{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploring sea ice heights with ICESat-2 (ATL07)\n",
    "\n",
    "Information obtained primarily from the ATL07/10 Algorithm Theoretical Basis Document (ATBD, Kwok et al., 2019) and the NSIDC product description page: https://nsidc.org/data/atl07.   \n",
    "\n",
    "* Notebook author: Alek Petty, relying extensively on above product manuals.    \n",
    "* Description: Notebook describing the ICESat-2 ATL07 product.   \n",
    "* Input requirements: Demo ATL07 data file   \n",
    "* Date: June 2019\n",
    "* More info: See the ATL07/ATL10 Algorithm Theoretical Basis Document (ATBD): https://icesat-2.gsfc.nasa.gov/sites/default/files/page_files/ICESat2_ATL07_ATL10_ATBD_r001.pdf    and the known issues document: https://nsidc.org/sites/nsidc.org/files/technical-references/ATL0710-KnownIssues.pdf\n",
    "\n",
    "\n",
    "## Notebook objectives\n",
    "* General understanding of the data included in a typical ATL07 file.\n",
    "* Reading in, plotting and basic analysis of ATL07 data.\n",
    "* How is ATL07 data used to generate ATL10 sea ice freeboards and what to look out for when using either product.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook instructions\n",
    "1. Follow along with the notebook tutorial. \n",
    "2. Play around changing options and re-running the relevant notebook cells. \n",
    "\n",
    "Here I use the HDF5 ATL07 file (ATL07-01_20181115003141_07240101_001_01.h5) from: https://nsidc.org/data/atl07. If using this using the ICESat-2 Pangeo instance, you can download the file from Amazon S3 using the notebook cell provided below.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ATL07 Background\n",
    "\n",
    "*NB: This is my laymans description of ATL07 compiled from the above resources, trying to condense this information to the key points of interest to potential sea ice users. let me know if you see any errors!*\n",
    "\n",
    "ATL07 is arguably the most important IS2 product for sea ice users. ATL07 provides along-track surface height and type (e.g. snow-covered ice, bare ice, open water) for the ice-covered seas of the northern and southern hemispheres. The primary input is the individual photon height estimates from the Level 2A Global Geolocated Photon Data (ATL03) product. Sea surface and sea ice height are estimated for segments along each of the six beams. \n",
    "\n",
    "The mean sea surface (MSS) derived from ICESat/CryoSat-2, the inverted barometric (IB) correction calculated from surface pressure from ATL09 and ocean tides are subtracted from the ATL03 surface heights (relative to WGS84). \n",
    "\n",
    "A two-step filtering process is then used to derive heights within a given segment (along-track series) of photons. First, A coarse surface filtering method is employed to remove obviously erroneuous returns from background/subsurface/clouds etc. Second, a fine surface filtering is applied which analyzes 150 photon segments to derive the mean height. A first-photon bias estimate is included in ATL07, based on system engineering with each height estimate. Subsurface-scattering, or volume scattering, is a bias that comes from photons that experience multiple scattering within the snow or ice before returning to the satellite.  This is not included in ATL07 but is something worth considering when using the data.\n",
    "\n",
    "Each of the 150 photon segments are then classified surface based on a decision-tree approach that determines the most likely surface type from three primary variables: the surface photon rate, the width of the photon distribution (or the fitted Gaussian) and the background rate. \n",
    "\n",
    "* The surface photon rate (photon returns per pulse) is a measure of the brightness, or apparent surface reflectance, of that height segment. In general, low surface rates indicate water or thin ice in open leads. \n",
    "* The width of the photon distribution provides a measure of the surface roughness and can be used to partition the height segments within four ranges that correspond to different surface types.\n",
    "* The background rate provides useful additional information when the solar elevation is high and sufficient photons are present. For example, a lack of relationship between sun angle and background photon rate  can indicates shadows (cloud shadows or ridge shadows), specular returns, or possibly, atmospheric effects. \n",
    "\n",
    "Much of this approach was based on data collected by NASA's MABEL (Multiple Altimeter Beam Experimental Lidar) photon-counting lidar prior to the ICESat-2 launch. Calibration/validation of the ATL07 data product is underway (Operation IceBridge spring 2019 data shortly to be made available).\n",
    "\n",
    "Below is an example of how one can read in and explore a given ATL07 file. Very open to suggestions on possible improvements!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Magic function to enable interactive plotting in Jupyter notebook\n",
    "#Allows you to zoom/pan within plots after generating\n",
    "#Normally, this would be %matplotlib notebook, but since we're using Juptyerlab, we need a different widget\n",
    "#%matplotlib notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "#Import necesary modules\n",
    "#Use shorter names (np, pd, plt) instead of full (numpy, pandas, matplotlib.pylot) for convenience\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import pandas as pd\n",
    "import numpy.ma as ma\n",
    "import h5py\n",
    "import s3fs\n",
    "import readers as rd\n",
    "import utils as ut\n",
    "import xarray as xr\n",
    "from astropy.time import Time\n",
    "\n",
    "# Use seasborn for nicer looking inline plots if available \n",
    "#import seaborn as sns\n",
    "#sns.set(context='notebook', style='darkgrid')\n",
    "#st = axes_style(\"whitegrid\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Beam selection   \n",
    "There are 6 beams to choose from in the ICESat-2 products (3 pairs of a strong and weak beam). The energy ratio between the weak and strong beams are  approximately 1:4 and are separated by 90 m in the across-track direction. The beam pairs are separated by ~3.3 km in the across-track direction, and the strong and weak beams are separated by ~2.5 km in the along-track direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "beamNum=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If file stored locally...\n",
    "#file_path = '../Data/'\n",
    "#ATL07_filename = 'ATL07-01_20181115003141_07240101_001_01.h5'\n",
    "#localPath = file_path + ATL07_filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If running on Pangeo instance and grabbing data from Amazon S3\n",
    "# Comment out the last command if you've already got the data in the Data dir\n",
    "bucket = 'pangeo-data-upload-oregon'\n",
    "fs = s3fs.S3FileSystem()\n",
    "dataDir = 'pangeo-data-upload-oregon/icesat2/'\n",
    "s3List = fs.ls(dataDir)\n",
    "#print(s3List)\n",
    "ATL07file='ATL07-01_20181115003141_07240101_001_01.h5'\n",
    "s3File='pangeo-data-upload-oregon/icesat2/'+ATL07file\n",
    "localFilePath='../Data/'+ATL07file\n",
    "#fs.get(s3File, localFilePath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getATL07data(fileT, numpy=False, beamNum=1, maxElev=1e6):\n",
    "    \"\"\" Pandas/numpy ATL07 reader\n",
    "    Written by Alek Petty, June 2018 (alek.a.petty@nasa.gov)\n",
    "\n",
    "    I've picked out the variables from ATL07 I think are of most interest to sea ice users, but by no means is this an exhastive list. \n",
    "    See the xarray or dictionary readers to load in the more complete ATL07 dataset\n",
    "    or explore the hdf5 files themselves (I like using the app Panpoly for this) to see what else you might want\n",
    "    \n",
    "    Args:\n",
    "        fileT (str): File path of the ATL07 dataset\n",
    "        numpy (flag): Binary flag for outputting numpy arrays (True) or pandas dataframe (False)\n",
    "        beamNum (int): ICESat-2 beam number (1 to 6)\n",
    "        maxElev (float): maximum surface elevation to remove anomalies\n",
    "\n",
    "    returns:\n",
    "        either: select numpy arrays or a pandas dataframe\n",
    "        \n",
    "    Updates:\n",
    "        V3 (June 2018) added observatory orientation flag, read in the beam number, not the string\n",
    "        V2 (June 2018) used astropy to more simply generate a datetime instance form the gps time\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    # Open the file\n",
    "    try:\n",
    "        ATL07 = h5py.File(fileT, 'r')\n",
    "    except:\n",
    "        return 'Not a valid file'\n",
    "    \n",
    "    #flag_values: 0, 1, 2; flag_meanings : backward forward transition\n",
    "    orientation_flag=ATL07['orbit_info']['sc_orient'][:]\n",
    "    \n",
    "    if (orientation_flag==0):\n",
    "        print('Backward orientation')\n",
    "        beamStrs=['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']\n",
    "                \n",
    "    elif (orientation_flag==1):\n",
    "        print('Forward orientation')\n",
    "        beamStrs=['gt3r', 'gt3l', 'gt2r', 'gt2l', 'gt1r', 'gt1l']\n",
    "        \n",
    "    elif (orientation_flag==2):\n",
    "        print('Transitioning, do not use for science!')\n",
    "    \n",
    "    beamStr=beamStrs[beamNum-1]\n",
    "    print(beamStr)\n",
    "    \n",
    "    lons=ATL07[beamStr+'/sea_ice_segments/longitude'][:]\n",
    "    lats=ATL07[beamStr+'/sea_ice_segments/latitude'][:]\n",
    "    \n",
    "    # Along track distance \n",
    "    # I removed the first point so it's distance relative to the start of the beam\n",
    "    along_track_distance=ATL07[beamStr+'/sea_ice_segments/seg_dist_x'][:] - ATL07[beamStr+'/sea_ice_segments/seg_dist_x'][0]\n",
    "    # Height segment ID (10 km segments)\n",
    "    height_segment_id=ATL07[beamStr+'/sea_ice_segments/height_segment_id'][:] \n",
    "    # Number of seconds since the GPS epoch on midnight Jan. 6, 1980 \n",
    "    delta_time=ATL07[beamStr+'/sea_ice_segments/delta_time'][:] \n",
    "    # Add this value to delta time parameters to compute full gps time\n",
    "    atlas_epoch=ATL07['/ancillary_data/atlas_sdp_gps_epoch'][:] \n",
    "\n",
    "    leapSecondsOffset=37\n",
    "    gps_seconds = atlas_epoch[0] + delta_time - leapSecondsOffset\n",
    "    # Use astropy to convert from gps time to datetime\n",
    "    tgps = Time(gps_seconds, format='gps')\n",
    "    tiso = Time(tgps, format='datetime')\n",
    "    \n",
    "    # Primary variables of interest\n",
    "    \n",
    "    # Beam segment height\n",
    "    elev=ATL07[beamStr+'/sea_ice_segments/heights/height_segment_height'][:]\n",
    "    # Flag for potential leads, 0=sea ice, 1 = sea surface\n",
    "    ssh_flag=ATL07[beamStr+'/sea_ice_segments/heights/height_segment_ssh_flag'][:] \n",
    "    \n",
    "    #Quality metrics for each segment include confidence level in the surface height estimate, \n",
    "    # which is based on the number of photons, the background noise rate, and the error measure provided by the surface-finding algorithm.\n",
    "    # Height quality flag, 1 for good fit, 0 for bad\n",
    "    quality=ATL07[beamStr+'/sea_ice_segments/heights/height_segment_quality'][:] \n",
    "    \n",
    "    elev_rms = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_rms'][:] #RMS difference between modeled and observed photon height distribution\n",
    "    seg_length = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_length_seg'][:] # Along track length of segment\n",
    "    height_confidence = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_confidence'][:] # Height segment confidence flag\n",
    "    reflectance = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_asr_calc'][:] # Apparent surface reflectance\n",
    "    ssh_flag = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_ssh_flag'][:] # Flag for potential leads, 0=sea ice, 1 = sea surface\n",
    "    seg_type = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_type'][:] # 0 = Cloud covered\n",
    "    gauss_width = ATL07[beamStr+'/sea_ice_segments/heights/height_segment_w_gaussian'][:] # Width of Gaussian fit\n",
    "\n",
    "    # Geophysical corrections\n",
    "    # NOTE: All of these corrections except ocean tides, DAC, \n",
    "    # and geoid undulations were applied to the ATL03 photon heights.\n",
    "    \n",
    "    # AVISO dynamic Atmospheric Correction (DAC) including inverted barometer (IB) effect (±5cm)\n",
    "    dac = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_dac'][:] \n",
    "    # Solid Earth Tides (±40 cm, max)\n",
    "    earth = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_earth'][:]\n",
    "    # Geoid (-105 to +90 m, max)\n",
    "    geoid = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_geoid'][:] \n",
    "    # Local displacement due to Ocean Loading (-6 to 0 cm)\n",
    "    loadTide = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_load'][:] \n",
    "    # Ocean Tides including diurnal and semi-diurnal (harmonic analysis), \n",
    "    # and longer period tides (dynamic and self-consistent equilibrium) (±5 m)\n",
    "    oceanTide = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_ocean'][:]\n",
    "    # Deformation due to centrifugal effect from small variations in polar motion \n",
    "    # (Solid Earth Pole Tide) (±1.5 cm, the ocean pole tide ±2mm amplitude is considered negligible)\n",
    "    poleTide = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_pole'][:] \n",
    "    # Mean sea surface (±2 m)\n",
    "    # Taken from ICESat and CryoSat-2, see Kwok and Morison [2015])\n",
    "    mss = ATL07[beamStr+'/sea_ice_segments/geophysical/height_segment_mss'][:]\n",
    "    \n",
    "    # Photon rate of the given segment\n",
    "    photon_rate = ATL07[beamStr+'/sea_ice_segments/stats/photon_rate'][:]\n",
    "    \n",
    "    # Estimated background rate from sun angle, reflectance, surface slope\n",
    "    background_rate = ATL07[beamStr+'/sea_ice_segments/stats/backgr_calc'][:]\n",
    "    \n",
    "    \n",
    "    \n",
    "    ATL07.close()\n",
    "    \n",
    "    if numpy:\n",
    "        # list the variables you want to output here..\n",
    "        return along_track_dist, elev\n",
    "    \n",
    "    else:\n",
    "        dF = pd.DataFrame({'elev':elev, 'lons':lons, 'lats':lats, 'ssh_flag':ssh_flag,\n",
    "                          'quality_flag':quality,\n",
    "                           'delta_time':delta_time,\n",
    "                           'along_track_distance':along_track_distance,\n",
    "                           'height_segment_id':height_segment_id, \n",
    "                           'photon_rate':photon_rate,'background_rate':background_rate,\n",
    "                          'datetime':tiso, 'mss': mss, 'seg_length':seg_length})\n",
    "        \n",
    "         # Add the datetime string\n",
    "        #dFtimepd=pd.to_datetime(dFtime)\n",
    "        #dF['datetime'] = pd.Series(dFtimepd, index=dF.index)\n",
    "        \n",
    "        # Filter out high elevation values \n",
    "        dF = dF[(dF['elev']<maxElev)]\n",
    "        # Reset row indexing\n",
    "        dF=dF.reset_index(drop=True)\n",
    "        return dF\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read in the data using the pandas reader above. Copied from the readers.py script.\n",
    "\n",
    "Take a look at the top few rows (change the number in head to increase this..)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward orientation\n",
      "gt3r\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>elev</th>\n",
       "      <th>lons</th>\n",
       "      <th>lats</th>\n",
       "      <th>ssh_flag</th>\n",
       "      <th>quality_flag</th>\n",
       "      <th>delta_time</th>\n",
       "      <th>along_track_distance</th>\n",
       "      <th>height_segment_id</th>\n",
       "      <th>photon_rate</th>\n",
       "      <th>background_rate</th>\n",
       "      <th>datetime</th>\n",
       "      <th>mss</th>\n",
       "      <th>seg_length</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.151697</td>\n",
       "      <td>-168.187226</td>\n",
       "      <td>73.235994</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2.747825e+07</td>\n",
       "      <td>2568.670525</td>\n",
       "      <td>10</td>\n",
       "      <td>0.780220</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2018-11-15 00:50:49.971690</td>\n",
       "      <td>-0.149128</td>\n",
       "      <td>121.853638</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.153673</td>\n",
       "      <td>-168.187239</td>\n",
       "      <td>73.236023</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2.747825e+07</td>\n",
       "      <td>2571.852364</td>\n",
       "      <td>11</td>\n",
       "      <td>1.516129</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2018-11-15 00:50:49.972142</td>\n",
       "      <td>-0.149039</td>\n",
       "      <td>59.295998</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.160038</td>\n",
       "      <td>-168.187248</td>\n",
       "      <td>73.236042</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.747825e+07</td>\n",
       "      <td>2574.083009</td>\n",
       "      <td>12</td>\n",
       "      <td>2.854167</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2018-11-15 00:50:49.972459</td>\n",
       "      <td>-0.148990</td>\n",
       "      <td>33.133297</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.167672</td>\n",
       "      <td>-168.187258</td>\n",
       "      <td>73.236063</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.747825e+07</td>\n",
       "      <td>2576.349800</td>\n",
       "      <td>13</td>\n",
       "      <td>5.071429</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2018-11-15 00:50:49.972780</td>\n",
       "      <td>-0.148929</td>\n",
       "      <td>19.062412</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.059221</td>\n",
       "      <td>-168.187298</td>\n",
       "      <td>73.236147</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2.747825e+07</td>\n",
       "      <td>2585.810232</td>\n",
       "      <td>14</td>\n",
       "      <td>2.607843</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2018-11-15 00:50:49.974118</td>\n",
       "      <td>-0.148674</td>\n",
       "      <td>35.393513</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       elev        lons       lats  ssh_flag  quality_flag    delta_time  \\\n",
       "0 -0.151697 -168.187226  73.235994         1             0  2.747825e+07   \n",
       "1 -0.153673 -168.187239  73.236023         1             0  2.747825e+07   \n",
       "2 -0.160038 -168.187248  73.236042         0             0  2.747825e+07   \n",
       "3 -0.167672 -168.187258  73.236063         0             0  2.747825e+07   \n",
       "4 -0.059221 -168.187298  73.236147         0             0  2.747825e+07   \n",
       "\n",
       "   along_track_distance  height_segment_id  photon_rate  background_rate  \\\n",
       "0           2568.670525                 10     0.780220              0.0   \n",
       "1           2571.852364                 11     1.516129              0.0   \n",
       "2           2574.083009                 12     2.854167              0.0   \n",
       "3           2576.349800                 13     5.071429              0.0   \n",
       "4           2585.810232                 14     2.607843              0.0   \n",
       "\n",
       "                     datetime       mss  seg_length  \n",
       "0  2018-11-15 00:50:49.971690 -0.149128  121.853638  \n",
       "1  2018-11-15 00:50:49.972142 -0.149039   59.295998  \n",
       "2  2018-11-15 00:50:49.972459 -0.148990   33.133297  \n",
       "3  2018-11-15 00:50:49.972780 -0.148929   19.062412  \n",
       "4  2018-11-15 00:50:49.974118 -0.148674   35.393513  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dF07=getATL07data(localFilePath, beamNum=beamNum)\n",
    "dF07.head(5)\n",
    "\n",
    "# Or get data using xarray\n",
    "#dF07X= rd.getATL07xarray(localFilePath, beamStr)\n",
    "#dF07X\n",
    "\n",
    "# ...Or data using numpy\n",
    "#along_track_dist, elev=getATL07data(ATL07_file_path, numpy=1, beam=beamStr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Map the data for visual inspection using Cartopy \n",
    "*NB (Basemap is often used for mapping but is not being officially supported by the community anymore)*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select variable of interest from the dataframe columns\n",
    "var='mss'\n",
    "dF07.along_track_distance\n",
    "plt.figure(figsize=(7,7), dpi= 90)\n",
    "# Make a new projection \"NorthPolarStereo\"\n",
    "ax = plt.axes(projection=ccrs.NorthPolarStereo(true_scale_latitude=70))\n",
    "plt.scatter(dF07['lons'], dF07['lats'],c=dF07[var], cmap='viridis', transform=ccrs.PlateCarree())\n",
    "#plt.pcolormesh(lons, lats, tile_to_plot,\n",
    "#               transform=ccrs.PlateCarree());\n",
    "\n",
    "ax.coastlines()\n",
    "#ax.drawmeridians()\n",
    "plt.colorbar(label=var, shrink=0.5, extend='both')\n",
    "\n",
    "# Limit the map to -60 degrees latitude and below.\n",
    "ax.set_extent([-180, 180, 90, 60], ccrs.PlateCarree())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get NESOSIM data which has variables including daily ice concentration from CDR dataset...\n",
    "\n",
    "dateStr='20181115'\n",
    "NESOSIMfilePath='NESOSIM-OSISAFsig150_ERAI_sf_SICCDR_Rhovariable_IC3_DYN1_WP1_LL1_WPF5.8e-07_WPT5_LLF2.9e-07-100kmnrt3-15082018-31012019.nc'\n",
    "localFilePath='../Data/'\n",
    "fs = s3fs.S3FileSystem()\n",
    "dataDir = 'pangeo-data-upload-oregon/icesat2/'\n",
    "#fs.get(dataDir+NESOSIMfilePath, localFilePath+NESOSIMfilePath)\n",
    "dNday= ut.getNESOSIM(localFilePath+NESOSIMfilePath, dateStr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Plot with ice concentration (or change to snowDepth) as a background\n",
    "\n",
    "# Select variable of interest from the dataframe columns\n",
    "var='photon_rate'\n",
    "dF07.along_track_distance\n",
    "plt.figure(figsize=(7,7), dpi= 90)\n",
    "# Make a new projection \"NorthPolarStereo\"\n",
    "ax = plt.axes(projection=ccrs.NorthPolarStereo(true_scale_latitude=70))\n",
    "plt.pcolormesh(dNday['longitude'], dNday['latitude'],ma.masked_where(dNday['iceConc']<0.3, dNday['iceConc']) , cmap='Blues_r', transform=ccrs.PlateCarree())\n",
    "\n",
    "plt.scatter(dF07['lons'], dF07['lats'],c=dF07[var], cmap='viridis', transform=ccrs.PlateCarree())\n",
    "#plt.pcolormesh(lons, lats, tile_to_plot,\n",
    "#               transform=ccrs.PlateCarree());\n",
    "\n",
    "ax.coastlines()\n",
    "#ax.drawmeridians()\n",
    "plt.colorbar(label=var, shrink=0.5, extend='both')\n",
    "\n",
    "# Limit the map to -60 degrees latitude and below.\n",
    "ax.set_extent([-180, 180, 90, 60], ccrs.PlateCarree())\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the segment heights of this section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot((dF07.along_track_distance)/1000., dF07.elev, color='r', marker='.', linestyle='None', alpha=0.2)\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Elevation (m)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Use the Pandas Groupby function to group the dataframe based on a given condition \n",
    "*surface type classification in this example..*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dFstype=dF07.groupby('ssh_flag')\n",
    "dFstype['elev'].agg(['mean', 'std', 'median', 'mad'])\n",
    "# Ice surface photons\n",
    "dFstypeIce=dFstype.get_group(0)\n",
    "# Lead/Sea surface photons\n",
    "dFstypeLeads=dFstype.get_group(1)\n",
    "\n",
    "# Note that in the ATL03 example I don't bother doing this \n",
    "# and just keep all the data in the table and just use a condition to display data \n",
    "# like the example commented out in the next cell..\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the grouped/classified data\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot((dFstypeIce.along_track_distance)/1000., dFstypeIce.elev, color='m', marker='.', linestyle='None', label='Sea ice', alpha=0.2)\n",
    "plt.plot((dFstypeLeads.along_track_distance)/1000., dFstypeLeads.elev, color='k', marker='x', linestyle='None',label='Leads', alpha=1.)\n",
    "plt.legend(frameon=False)\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Elevation (m)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An alternative (more intuitive?) approach...\n",
    "\n",
    "#plt.figure(figsize=(12, 5))\n",
    "#plt.plot((dF07[(dF07['ssh_flag']>0)]['along_track_distance']-dF07['along_track_distance'][0])/1000., dF07[(dF07['ssh_flag']>0)]['elev'], color='k', marker='x', linestyle='None',label='sea surface', alpha=1)\n",
    "#plt.plot((dF07[(dF07['ssh_flag']==0)]['along_track_distance']-dF07['along_track_distance'][0])/1000., dF07[(dF07['ssh_flag']==0)]['elev'], color='b', marker='.', linestyle='None',label='ice', alpha=0.3)\n",
    "#plt.legend(frameon=False)\n",
    "#plt.xlabel('Along track distance (km)')\n",
    "#plt.ylabel('Elevation (m)')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's simplify things by just looking at a 10 km along-track section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Section of ATL07 data\n",
    "sectionNum=200 # (NB 220 is a nice example, 200 seems to show surprisngly bad lead precisions..?)\n",
    "# Section length (for plotting purposes) in meters\n",
    "sectionSize=10000.\n",
    "\n",
    "# Find data that satisfies these conditions and then group the data like before\n",
    "idx=np.where((dF07['along_track_distance']>sectionNum*sectionSize)&(dF07['along_track_distance']<(sectionNum+1)*10000.))[0]\n",
    "df07S=dF07.iloc[idx]\n",
    "dFStype=df07S.groupby('ssh_flag')\n",
    "dFStype['elev'].agg(['mean', 'std', 'median', 'mad'])\n",
    "dFStypeIce=dFStype.get_group(0)\n",
    "dFStypeLeads=dFStype.get_group(1)\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot((dFStypeIce.along_track_distance)/1000., dFStypeIce.elev, color='m', marker='.', linestyle='None', label='Sea ice', alpha=0.2)\n",
    "plt.plot((dFStypeLeads.along_track_distance)/1000., dFStypeLeads.elev, color='k', marker='x', linestyle='None',label='Leads', alpha=1.)\n",
    "plt.legend(frameon=False)\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Elevation (m)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's calculate the mean elevation of the sea surface (and sea ice) as a simple test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sea ice elevation\n",
    "meanIceElev=dFstype['elev'].get_group(0).mean()\n",
    "# Sea surface elevation\n",
    "meanSSH=dFstype['elev'].get_group(1).mean()\n",
    "print('Sea ice elevation (m):', meanIceElev)\n",
    "print('SSH (m):', meanSSH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### OK well now it's clearly pretty simple to derive some freeboard!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot((dFStypeIce.along_track_distance)/1000., dFStypeIce.elev-meanSSH, color='g', marker='.', linestyle='-', label='Freeboard', alpha=0.5)\n",
    "#plt.legend(frameon=False)\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Freeboard (m)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explore the photon rate\n",
    "*Note the higher photon rates where we have leads!*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 5))\n",
    "plt.scatter((dFStypeIce.along_track_distance)/1000., dFStypeIce.elev-meanSSH, c=dFStypeIce.photon_rate, label='photon_rate')\n",
    "#plt.legend(frameon=False)\n",
    "plt.colorbar(label='Photon rate')\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Freeboard (m)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explore the background rate\n",
    "*Note that this is the calculated background rate in ATL07 based on the sun angle, surface slope, unit reflectance. There is also an observed background rate at a given along-track posting (25 Hz or 200 Hz)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot((dFStypeIce.along_track_distance)/1000., dFStypeIce.background_rate, color='k', marker='.', linestyle='-', label='background_rate', alpha=0.5)\n",
    "#plt.legend(frameon=False)\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Background rate (Hz)')\n",
    "plt.show()\n",
    "\n",
    "# Looks like it's nighttime in the Arctic!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 5))\n",
    "plt.plot((dFStypeIce.along_track_distance)/1000., dFStypeIce.seg_length, color='g', marker='.', linestyle='', label='segment length', alpha=0.5)\n",
    "#plt.legend(frameon=False)\n",
    "plt.xlabel('Along track distance (km)')\n",
    "plt.ylabel('Segment length (m)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extra ideas\n",
    "\n",
    "1. Try downloading some more ATL07 data from the NSIDC (following the hackweek tutorial) and see what it looks like when using it in this processing chain. \n",
    "2. Explore the photon classification scheme.\n",
    "3. Explore the photon rate and background rate. How do they variable with the open water/sea ice classification? Do some scatter plots of photon rate versus ice type. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Onwards to the ATL10 Notebook...!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
